    module math_mod
    use spark_cubed_sphere, only: i_kind, r_kind, r8
    contains
    subroutine calc_polynomial_matrix(d,m,n,xi,eta,A)
      integer(i_kind), intent(in ) :: d ! polynomial degree
      integer(i_kind), intent(in ) :: m ! number of points
      integer(i_kind), intent(in ) :: n ! number of terms of polynomial, n = (d+1)*(d+2)/2
      real   (r_kind), intent(in ) :: xi (m)
      real   (r_kind), intent(in ) :: eta(m)
      real   (r_kind), intent(out) :: A  (m,n)
      
      real   (r_kind) :: x
      real   (r_kind) :: y
      integer(i_kind) :: iPOC
      integer(i_kind)  :: i,j,k
      
      do iPOC = 1,m
        x = xi (iPOC)
        y = eta(iPOC)
        
        k = 0
        do j = 0,d
          do i = 0,j
            k = k + 1
            A(iPOC,k) = x**real(j-i,r_kind)*y**real(i,r_kind)
          enddo
        enddo
      enddo
  
    end subroutine calc_polynomial_matrix
    
    subroutine calc_polynomial_deriv_matrix(d,m,n,xi,eta,px,py)
      integer(i_kind), intent(in ) :: d ! polynomial degree
      integer(i_kind), intent(in ) :: m ! number of points
      integer(i_kind), intent(in ) :: n ! number of points on cell, n = (d+1)*(d+2)/2
      real   (r_kind), intent(in ) :: xi (m)
      real   (r_kind), intent(in ) :: eta(m)
      real   (r_kind), intent(out) :: px (m,n)
      real   (r_kind), intent(out) :: py (m,n)
      
      real(16) x
      real(16) y
      integer  iPOC, i, j, k
      real(16) powx1,powx2
      real(16) powy1,powy2
      real(16) coefx,coefy
      
      do iPOC = 1,m
        x = xi (iPOC)
        y = eta(iPOC)
        
        k = 0
        do j = 0,d
          do i = 0,j
            k = k + 1
            powx1 = merge( 0._r_kind, real(j-i-1,r_kind), real(j-i-1,r_kind)<0._r_kind )
            powx2 = real(i  ,r_kind)
            powy1 = real(j-i,r_kind)
            powy2 = merge( 0._r_kind, real(i-1  ,r_kind), real(i-1  ,r_kind)<0._r_kind )
            coefx = merge( 0._r_kind, real(j-i  ,r_kind), real(j-i-1,r_kind)<0._r_kind )
            coefy = merge( 0._r_kind, real(i    ,r_kind), real(i-1  ,r_kind)<0._r_kind )
            
            px(iPOC,k) = coefx * x**powx1 * y**powx2
            py(iPOC,k) = coefy * x**powy1 * y**powy2
          enddo
        enddo
      enddo
    
    end subroutine calc_polynomial_deriv_matrix
  
    subroutine  calc_polynomial_triangle_integration(d,c)
      integer(i_kind), intent(in ) :: d ! degree of polynomial
      real   (r_kind), intent(out) :: c(:)
      
      integer(i_kind) :: i,j,k,r
      
      k = 0
      c = 0
      do j = 0,d
        do i = 0,j
          k = k + 1
          do r = 0,i+1
            c(k) = c(k) + (-1)**(r) * nchoosek(i+1,r) / real( j - i + r + 1, r_kind )
          enddo
          c(k) = c(k) / real( i + 1, r_kind )
        enddo
      enddo
      
    end subroutine  calc_polynomial_triangle_integration
    
    subroutine  calc_polynomial_square_integration(d,x_min,x_max,y_min,y_max,c)
      integer(i_kind), intent(in ) :: d ! degree of polynomial
      real   (r_kind), intent(in ) :: x_min
      real   (r_kind), intent(in ) :: x_max
      real   (r_kind), intent(in ) :: y_min
      real   (r_kind), intent(in ) :: y_max
      real   (r_kind), intent(out) :: c(:)
      
      integer(i_kind) :: i,j,k
      
      k = 0
      c = 0
      do j = 0,d
        do i = 0,j
          k = k + 1
          c(k) = ( x_max**(j-i+1) - x_min**(j-i+1) ) * ( y_max**(i+1) - y_min**(i+1) ) / real( ( i + 1 ) * ( j - i + 1 ), r_kind )
        enddo
      enddo
      
    end subroutine  calc_polynomial_square_integration
    
    subroutine calc_rectangle_poly_matrix(nx,ny,m,xi,eta,A,existPolyTerm)
      integer(i_kind), intent(in   ) :: nx ! number of points on x direction for reconstruction
      integer(i_kind), intent(in   ) :: ny ! number of points on y direction for reconstruction
      integer(i_kind), intent(in   ) :: m  ! number of unkonwn point values
      real   (r_kind), intent(in   ) :: xi (m)
      real   (r_kind), intent(in   ) :: eta(m)
      real   (r_kind), intent(inout) :: A  (:,:)
      real   (r_kind), intent(in   ),optional :: existPolyTerm(nx*ny)
      
      real   (r_kind) :: ext(nx*ny)
      
      real   (r_kind) :: x
      real   (r_kind) :: y
      integer(i_kind) :: iPOC
      integer(i_kind)  :: i,j,k,iCOS
      
      ext = 1
      if(present(existPolyTerm))ext = existPolyTerm
      
      do iPOC = 1,m
        x = xi (iPOC)
        y = eta(iPOC)
        
        k = 0
        iCOS = 0
        do j = 0,ny-1
          do i = 0,nx-1
            k = k + 1
            if(ext(k)>0)then
              iCOS = iCOS + 1
              A(iPOC,iCOS) = x**real(i,r_kind) * y**real(j,r_kind)
            endif
          enddo
        enddo
      enddo
  
    end subroutine calc_rectangle_poly_matrix
    
    subroutine calc_rectangle_poly_integration(nx,ny,x_min,x_max,y_min,y_max,c,existPolyTerm)
      integer(i_kind), intent(in   ) :: nx  ! number of points on x direction
      integer(i_kind), intent(in   ) :: ny  ! number of points on y direction
      real   (r_kind), intent(in   ) :: x_min
      real   (r_kind), intent(in   ) :: x_max
      real   (r_kind), intent(in   ) :: y_min
      real   (r_kind), intent(in   ) :: y_max
      real   (r_kind), intent(inout) :: c(:)
      real   (r_kind), intent(in   ),optional :: existPolyTerm(nx*ny)
      
      real   (r_kind) :: ext(nx*ny)
      integer(i_kind) :: i,j,k,iCOS
      
      ext = 1
      if(present(existPolyTerm))ext = existPolyTerm
              
      k    = 0
      c    = 0
      iCOS = 0
      do j = 0,ny-1
        do i = 0,nx-1
          k = k + 1
          if(ext(k)>0)then
            iCOS = iCOS + 1
            c(iCOS) = ( x_max**(i+1) - x_min**(i+1) ) * ( y_max**(j+1) - y_min**(j+1) ) / real( ( i + 1 ) * ( j + 1 ), r_kind )
          endif
        enddo
      enddo
      
    end subroutine  calc_rectangle_poly_integration
    
    subroutine calc_rectangle_poly_deriv_matrix(nx,ny,m,xi,eta,dpdx,dpdy,existPolyTerm)
      integer(i_kind), intent(in   ) :: nx ! number of points on x direction for reconstruction
      integer(i_kind), intent(in   ) :: ny ! number of points on y direction for reconstruction
      integer(i_kind), intent(in   ) :: m  ! number of unkonwn point values
      real   (r_kind), intent(in   ) :: xi (m)
      real   (r_kind), intent(in   ) :: eta(m)
      real   (r_kind), intent(inout) :: dpdx(:,:)
      real   (r_kind), intent(inout) :: dpdy(:,:)
      real   (r_kind), intent(in   ),optional :: existPolyTerm(nx*ny)
      
      real   (r_kind) :: ext(nx*ny)
      
      real   (r_kind) :: x
      real   (r_kind) :: y
      integer(i_kind) :: iPOC
      integer(i_kind)  :: i,j,k,iCOS
      
      ext = 1
      if(present(existPolyTerm))ext = existPolyTerm
      
      dpdx = 0
      dpdy = 0
      do iPOC = 1,m
        x = xi (iPOC)
        y = eta(iPOC)
        
        k = 0
        iCOS = 0
        do j = 0,ny-1
          do i = 0,nx-1
            k = k + 1
            if(ext(k)>0)then
              iCOS = iCOS + 1
              dpdx(iPOC,iCOS) = merge(0.0_r_kind, real(i,r_kind),i-1<0) * x**merge(0,i-1,i-1<0) * y**j
              dpdy(iPOC,iCOS) = merge(0.0_r_kind, real(j,r_kind),j-1<0) * x**i * y**merge(0,j-1,j-1<0)
            endif
          enddo
        enddo
      enddo
  
    end subroutine calc_rectangle_poly_deriv_matrix
    
    ! calculate nth derivative for rectangular polynomial
    subroutine calc_rectangle_poly_deriv_n_matrix(nx,ny,m,xi,eta,dp,existPolyTerm)
      integer(i_kind), intent(in   ) :: nx ! number of points on x direction for reconstruction
      integer(i_kind), intent(in   ) :: ny ! number of points on y direction for reconstruction
      integer(i_kind), intent(in   ) :: m  ! number of unkonwn point values
      real   (r_kind), intent(in   ) :: xi (m)
      real   (r_kind), intent(in   ) :: eta(m)
      real   (r_kind), intent(inout) :: dp(0:nx-1,0:ny-1,m,nx*ny)
      real   (r_kind), intent(in   ),optional :: existPolyTerm(nx*ny)
      
      real   (r_kind) :: ext(nx*ny)
      
      real   (r_kind) :: x
      real   (r_kind) :: y
      integer(i_kind) :: iPOC
      integer(i_kind) :: id,jd
      integer(i_kind) :: i,j,k,iCOS
      
      real   (r_kind) :: cx(nx)
      real   (r_kind) :: cy(ny)
      
      do i = 1,nx
        cx(i) = i
      enddo
      do j = 1,ny
        cy(j) = j
      enddo
      
      ext = 1
      if(present(existPolyTerm))ext = existPolyTerm
      
      dp = 0
      do iPOC = 1,m
        x = xi (iPOC)
        y = eta(iPOC)
        
        k = 0
        iCOS = 0
        do j = 0,ny-1
          do i = 0,nx-1
            k = k + 1
            if(ext(k)>0)then
              iCOS = iCOS + 1
              do jd = 0,ny-1
                do id = 0,nx-1
                  dp(id,jd,iPOC,iCOS) = merge(1.0_r_kind, merge(0.0_r_kind, product(cx(nx-1:nx-id:-1)),i-id<0), id==0) * x**merge(0,i-id,i-id<0) &
                                      * merge(1.0_r_kind, merge(0.0_r_kind, product(cy(ny-1:ny-jd:-1)),j-jd<0), jd==0) * y**merge(0,j-jd,j-jd<0) &
                                      * existPolyTerm(k)
                enddo
              enddo
            endif
          enddo
        enddo
      enddo
  
    end subroutine calc_rectangle_poly_deriv_n_matrix
    
    ! spherical distance on unit sphere
    function spherical_distance(lat1,lon1,lat2,lon2,r)
      real(r_kind) :: spherical_distance
      real(r_kind),intent(in) :: lat1,lon1,lat2,lon2
      real(r_kind),intent(in) :: r
      
      !spherical_distance = r * acos( sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2)*cos(lon1-lon2) )
      spherical_distance = r * acos(min(1.0d0, max(-1.0d0, sin(lat1) * sin(lat2) + cos(lat1) * cos(lat2) * cos(lon1 - lon2))))
    end function spherical_distance
  
    function nchoosek(n,k) ! same as nchoosek in matlab
      real   (r_kind) :: nchoosek
      integer(i_kind) :: n
      integer(i_kind) :: k
      
      nchoosek = factorial(n) / ( factorial(n-k) * factorial(k) )
      
    end function nchoosek
    
    function factorial(n)
      real   (r_kind) :: factorial
      integer(i_kind) :: n
      
      factorial = gamma(real(n+1,r_kind))
    
    end function factorial
    
    ! calculate inverse matrix of A_input
    ! N is the order of matrix A_input and A
    ! A is inverse A_input
    ! L is status info
    SUBROUTINE BRINV(N,A_input,A,L)
    implicit none
    integer(i_kind),intent(in )           :: N
    real   (r_kind),intent(in )           :: A_input(N,N)
    real   (r_kind),intent(out)           :: A      (N,N)
    integer(i_kind),intent(out), optional :: L
    
    real   (r_kind) :: T,D
    integer(i_kind) :: IS(N),JS(N)
    integer(i_kind) :: i,j,k
    
    A = A_input
    
    if(present(L))L=1
    do K=1,N
      D=0._r_kind
      do I=K,N
        do J=K,N
          IF (ABS(A(I,J)).GT.D) THEN
            D=ABS(A(I,J))
            IS(K)=I
            JS(K)=J
          END IF
        enddo
      enddo
    
      IF (D+1._r_kind.EQ.1._r_kind) THEN
        if(present(L))L=0
        WRITE(*,*)'ERR**NOT INV'
        RETURN
      END IF
      
      do J=1,N
        T=A(K,J)
        A(K,J)=A(IS(K),J)
        A(IS(K),J)=T
      enddo
      
      do I=1,N
        T=A(I,K)
        A(I,K)=A(I,JS(K))
        A(I,JS(K))=T
      enddo
      
      A(K,K)=1._r_kind/A(K,K)
      do J=1,N
        IF (J.NE.K) THEN
          A(K,J)=A(K,J)*A(K,K)
        END IF
      enddo
      
      do I=1,N
        IF (I.NE.K) THEN
          do J=1,N
            IF (J.NE.K) THEN
              A(I,J)=A(I,J)-A(I,K)*A(K,J)
            END IF
          enddo
        END IF
      enddo
      
      do I=1,N
        IF (I.NE.K) THEN
          A(I,K)=-A(I,K)*A(K,K)
        END IF
      enddo
    enddo
    
    do K=N,1,-1
      do J=1,N
        T=A(K,J)
        A(K,J)=A(JS(K),J)
        A(JS(K),J)=T
      enddo
      do I=1,N
        T=A(I,K)
        A(I,K)=A(I,IS(K))
        A(I,IS(K))=T
      enddo
    enddo
    RETURN
    END SUBROUTINE BRINV
    
    pure real(r_kind) function center_diff2(f, dx) result(res)
    
      real(r_kind), intent(in) :: f(-1:1)
      real(r_kind), intent(in) :: dx
    
      res = (f(1) - f(-1)) / (2 * dx)
    
    end function center_diff2
    
    pure real(r_kind) function center_diff4(f, dx) result(res)
    
      real(r_kind), intent(in) :: f(-2:2)
      real(r_kind), intent(in) :: dx
    
      res = (f(-2) - 8 * f(-1) + 8 * f(1) - f(2)) / (12 * dx)
    
    end function center_diff4
    
    pure real(r_kind) function center_diff6(f, dx) result(res)
    
      real(r_kind), intent(in) :: f(-3:3)
      real(r_kind), intent(in) :: dx
    
      res = (-f(-3) + 9 * f(-2) - 45 * f(-1) + 45 * f(1) - 9 * f(2) + f(3)) / (60 * dx)
    
    end function center_diff6
    
    pure real(r_kind) function center_diff8(f, dx) result(res)
    
      real(r_kind), intent(in) :: f(-4:4)
      real(r_kind), intent(in) :: dx
      real(r_kind), dimension(9), parameter :: coef = (/1./280,-4./105,1./5,-4./5,0.,4./5,-1./5,4./105,-1./280/)

    
      res = dot_product(f,coef) / dx
    
    end function center_diff8
    
    pure real(r_kind) function center_diff12(f, dx) result(res)
    
      real(r_kind), intent(in) :: f(-6:6)
      real(r_kind), intent(in) :: dx
      real(r_kind), dimension(13), parameter :: coef = (/1./5544,-1./385,1./56,-5./63,15./56,-6./7,0.,6./7,-15./56,5./63,-1./56,1./385,-1./5544/)

    
      res = dot_product(f,coef) / dx
    
    end function center_diff12

  subroutine convert_hor_deriv_cube_to_sph(dzdlon,dzdlat,dzdx,dzdy,iA)
    real(r_kind), intent(out) :: dzdlon
    real(r_kind), intent(out) :: dzdlat
    real(r_kind), intent(in ) :: dzdx
    real(r_kind), intent(in ) :: dzdy
    real(r_kind), intent(in ) :: iA(2,2) ! inverse A matrix
    
    dzdlon = iA(1,1) * dzdx + iA(2,1) * dzdy
    dzdlat = iA(1,2) * dzdx + iA(2,2) * dzdy
    
  end subroutine convert_hor_deriv_cube_to_sph

  subroutine convert_hor_deriv_sph_to_cube(dzdx,dzdy,dzdlon,dzdlat,A)
    real(r_kind), intent(out) :: dzdx
    real(r_kind), intent(out) :: dzdy
    real(r_kind), intent(in ) :: dzdlon
    real(r_kind), intent(in ) :: dzdlat
    real(r_kind), intent(in ) :: A(2,2) ! A matrix
    
    dzdx = A(1,1) * dzdlon + A(2,1) * dzdlat
    dzdy = A(1,2) * dzdlon + A(2,2) * dzdlat
    
  end subroutine convert_hor_deriv_sph_to_cube
    
  subroutine interp_1d(sx,sy,tx,ty,dydx)
    real(r_kind), intent(in )          :: sx  (:) ! coordinate on source point
    real(r_kind), intent(in )          :: sy  (:) ! data value on source point
    real(r_kind), intent(in )          :: tx  (:) ! coordinate on target point
    real(r_kind), intent(out)          :: ty  (:) ! data value on target point
    real(r_kind), intent(out),optional :: dydx(:)
    
    integer :: ids,ide
    integer :: its,ite
    integer :: refIdx(1)
    integer :: refL, refR
    integer :: i
    logical :: flip_data = .false.
    
    real(r_kind) :: x(lbound(sx,1):ubound(sx,1))
    real(r_kind) :: y(lbound(sy,1):ubound(sy,1))
    real(r_kind) :: dx, coefL, coefR
    
    ids = lbound(x , 1)
    ide = ubound(x , 1)
    its = lbound(tx, 1)
    ite = ubound(tx, 1)
    
    x = sx
    y = sy
    
    flip_data = .false.
    if( x(2) - x(1) < 0 )then
      flip_data = .true.
      call flipud(x)
      call flipud(y)
    endif
    
    do i = its, ite
      ! Find the location index of target point
      refIdx = maxloc( x, x < tx(i) )
      
      if(refIdx(1)==0  ) refIdx(1) = 1       ! exterpolation
      if(refIdx(1)==ide) refIdx(1) = ide - 1 ! exterpolation
      
      refL = refIdx(1)
      refR = refL + 1
      
      dx = x(refR) - x(refL)
      
      coefL = (x (refR) - tx(i   )) / dx
      coefR = (tx(i   ) - x (refL)) / dx
      
      ty(i) = coefL * y(refL) + coefR * y(refR)
      
      if(present(dydx)) dydx(i) = (y(refR) - y(refL)) / dx
    enddo
    
  end subroutine interp_1d
  
  real(r_kind) function bilinear_interp(tx,ty,xL,xR,yB,yT,vBL,vBR,vTL,vTR,missing_value)
    real(r_kind), intent(in)          :: tx, ty          ! x and y coordinates of target point
    real(r_kind), intent(in)          :: xL,xR           ! x coordinate of source points, xL: left, xR: rigth
    real(r_kind), intent(in)          :: yB,yT           ! y coordinate of source points, yB: bottom, yT: top
    real(r_kind), intent(in)          :: vBL,vBR,vTL,vTR ! point values of source points, BL: bottom left, BR: bottom right, TL: up top, TR: top right
    real(r_kind), intent(in),optional :: missing_value   ! missing value
    
    real(r_kind) fB,fT
    real(r_kind) dx, dy
    real(r_kind) dL, dR, dB, dT
    real(r_kind) msv
    real(r_kind) max_value, min_value
    real(r_kind), parameter :: tol = 1.e-30
    
    if(present(missing_value))then
      msv = missing_value
    else
      msv = huge(msv)
    endif
    
    dx = xR - xL
    dy = yT - yB
    dL = tx - xL
    dR = xR - tx
    dB = ty - yB
    dT = yT - ty
    
    if(vBL/=msv.and.&
       vBR/=msv.and.&
       vTL/=msv.and.&
       vTR/=msv)then
      if( dx /= 0 )then
        fB = ( dR * vBL + dL* vBR ) / dx
        fT = ( dR * vTL + dL* vTR ) / dx
      else
        fB = vBL
        fT = vTL
      endif
      
      if( dy /=0 )then
        bilinear_interp = ( dB * fT + dT * fB ) / dy
      else
        bilinear_interp = fT
      endif
    else
      bilinear_interp = nearest_interp(tx,ty,xL,yB,xR,yB,xL,yT,xR,yT,vBL,vBR,vTL,vTR,missing_value)
    endif
    
    max_value = max(vBL,vBR,vTL,vTR)
    min_value = min(vBL,vBR,vTL,vTR)
    if( bilinear_interp > max_value .or. bilinear_interp < min_value )then
      if( bilinear_interp > max_value .and. max_value /= 0 .and. abs( ( bilinear_interp - max_value ) / max_value ) > 1.e-14 .and. abs(max_value) > tol  )then
        print*,'bilinear_interp, max_value, tx, ty, xL, xR, yB, yT',bilinear_interp, max_value, tx, ty, xL, xR, yB, yT
        stop 'Error in bilinear interpolation! output value is greater than sources'
      endif
      if( bilinear_interp < min_value .and. min_value /= 0 .and. abs( ( bilinear_interp - min_value ) / min_value ) > 1.e-14 .and. abs(min_value) > tol )then
        print*,'bilinear_interp, min_value, tx, ty, xL, xR, yB, yT',bilinear_interp, min_value, tx, ty, xL, xR, yB, yT
        stop 'Error in bilinear interpolation! output value is smaller than sources'
      endif
    endif
  
  end function bilinear_interp
    
  real(r_kind) function nearest_interp(lont,latt,lonBL,latBL,lonBR,latBR,lonTL,latTL,lonTR,latTR,vBL,vBR,vTL,vTR,missing_value)
    real(r_kind), intent(in)          :: lont ,latt              ! longitude and latitude of target point
    real(r_kind), intent(in)          :: lonBL,lonBR,lonTL,lonTR ! longitude of source points, BL: bottom left, BR: bottom right, TL: up top, TR: top right
    real(r_kind), intent(in)          :: latBL,latBR,latTL,latTR ! latitude of source points, BL: bottom left, BR: bottom right, TL: up top, TR: top right
    real(r_kind), intent(in)          :: vBL,vBR,vTL,vTR         ! point values of source points, BL: bottom left, BR: bottom right, TL: up top, TR: top right
    real(r_kind), intent(in),optional :: missing_value           ! missing value
    ! Attension, aBL of the longitudes and latitdues must be set in radian.
    
    real(r_kind) :: rBL,rBR,rTL,rTR   ! distance between target point and source points, BL: bottom left, BR: bottom right, TL: up top, TR: top right
    real(r_kind) :: radius = 6371229. ! radius of earth, just for computing the spherical distance, it does not influent the result of nearest interpolation
    real(r_kind) :: msv
    
    if(present(missing_value))then
      msv = missing_value
    else
      msv = huge(msv)
    endif
    
    rBL = radius * acos(sin(lont)*sin(lonBL)+cos(lont)*cos(lonBL)*cos(latt-latBL))
    rBR = radius * acos(sin(lont)*sin(lonBR)+cos(lont)*cos(lonBR)*cos(latt-latBR))
    rTL = radius * acos(sin(lont)*sin(lonTL)+cos(lont)*cos(lonTL)*cos(latt-latTL))
    rTR = radius * acos(sin(lont)*sin(lonTR)+cos(lont)*cos(lonTR)*cos(latt-latTR))
    
    if(            vBL/=msv) nearest_interp = vBL
    if(rBR<rBL.and.vBR/=msv) nearest_interp = vBR
    if(rTL<rBR.and.vTL/=msv) nearest_interp = vTL
    if(rTR<rTL.and.vTR/=msv) nearest_interp = vTR
  end function nearest_interp
  
  ! cross product for 3D vector
  function cross(a, b)
    real(r_kind), DIMENSION(3)             :: cross
    real(r_kind), DIMENSION(3), INTENT(IN) :: a, b
  
    cross(1) = a(2) * b(3) - a(3) * b(2)
    cross(2) = a(3) * b(1) - a(1) * b(3)
    cross(3) = a(1) * b(2) - a(2) * b(1)
  end function cross
  
  subroutine spline_integration(n,x,y,da,db,nt,t,ty)
    !---------------------------------subroutine  comment
    !
    !  Purpose   :   ������������֮��֧����
    !    
    !-----------------------------------------------------
    !  Input  parameters  :
    !       1.   n-----��ֵ�ڵ������1�����оŸ��ڵ� ��N=8
    !       2.   x ---�ڵ��Ա���  Ϊ��0��N��ά����
    !       3.   y----�ڵ������  ��0��N��ά����
    !       4.   nt Ҫ����������ά��
    !       5.   t Ҫ���������  (1:nt) ά����
    !       6.   da  ------��㴦���� f'(x(0))
    !       7.   db -------�յ㴦���� f'(x(N))
    !  Output parameters  :
    !       1.   ty ---Ҫ�������������1��nt��ά
    !
    !  Common parameters  :
    !
    !----------------------------------------------------
    !  Post Script :
    !            
    !            ����Ҫ��ֵ���������������Բ��ذ���С����
    !----------------------------------------------------
    
    implicit real(r_kind)(a-z)
    
    integer::n, nt
    !nΪ��ֵ�ڵ�����1���������9���ڵ㣬��n=8
    !nt Ҫ��ֵ��ĸ�����������t,ty��ά��
    
    integer::i,j,k
    
    real(r_kind)::x(0:n),y(0:n)
    
    real(r_kind)::t(nt),ty(nt)
    
    real(r_kind)::h(0:n-1)
    
    real(r_kind)::f1(0:n-1),f2(1:n-1)
    
    real(r_kind)::u(1:n-1),lambda(1:n-1),d(0:n)
    
    real(r_kind)::M(0:n)
    
    real(r_kind)::A(0:n,0:n)
    
    integer :: node(nt)  ! ��¼ÿ��������������ĸ��ڵ㷶Χ�ڣ���ڵ�1���ڵ�2֮�䶼�ǽڵ�1�Ŀ��Ʒ�Χ
    integer :: dnode(nt) ! ��¼ÿ�������֮����˶��ٸ��ڵ�
    
    do i=0,n-1
      h (i) = x(i+1) - x(i)
      f1(i) = ( y(i+1) - y(i) ) / h(i)
    end do
    
    !�Թ�֧�߽���������   ����  d(0) �� d(n)
    d(0) = 6./h(0  ) * (f1(0) - da     )
    d(n) = 6./h(n-1) * (db    - f1(n-1))
    
    !��� u, lambda, d
    do i=1,n-1
      u     (i) = h(i-1) / ( h(i-1) + h(i) )
      lambda(i) = 1 - u(i)
      
      f2(i) = ( f1(i-1) - f1(i) ) / ( x(i-1) - x(i+1) )
      
      d(i) = 6. * f2(i) 
    end do
    
    !����A����ֵ
    A = 0
    do i = 1, n-1
      a(i,i) = 2.
    end do
    
    do i = 2, n-1
      a(i,i-1) = u(i)
    end do
    
    do i = 1, n-2
      a(i,i+1) = lambda(i)
    end do
    
    !-------�������Ȼ������������Ҫ����A������ĩ��Ԫ�أ�������ͬ
    !���� A���������Ԫ��
    a(0,0) = 2.
    a(0,1) = 1.
    
    !����A����Ԫ�ص�ĩ��Ԫ��
    a(n,n-1) = 1.
    a(n,n  ) = 2.
    
    ! ����������ֵ
    d(1  ) = d(1  ) - u     (1  ) * M(0)
    d(n-1) = d(n-1) - lambda(n-1) * M(n)
    
    
    !call gauss(a,d,M,N+1)
    call chase(a,d,M,N+1)
    
    ! ���¿�ʼ�������ֵ
    do k=1,nt
      !------------
      !  ��Ҫ��ֵ����ÿ���������ԣ����ҵ����������е�λ��
      do i=1,n-1
        if (t(k)<x(i+1)) exit
      end do
      
      if(i==n)then
        ! ���
        i = n - 1
      endif
      
      node(k) = i
    enddo
    dnode(2:nt) = node(2:nt) - node(1:nt-1)
    dnode(1   ) = 0
    
    !�ֶλ���
    ty    = 0.
    k     = 1
    ty(k) = 0.
    do k = 2, nt
      if(dnode(k)==0)then
        ty(k) = ty(k) + intfunc( t(k),node(k) ) - intfunc( t(k-1),node(k) )
      else
        ty(k) = ty(k) + intfunc( x(node(k-1)+1),node(k-1) ) - intfunc( t(k-1),node(k-1) )
        
        if(dnode(k)>1)then
          do i = 1, dnode(k)
            ty(k) = ty(k) + intfunc( x(node(k-1)+i+1),node(k-1)+i ) - intfunc( x(node(k-1)+i),node(k-1)+i )
          enddo
        endif
        
        ty(k) = ty(k) + intfunc( t(k),node(k) ) - intfunc( x(node(k)),node(k) )
      endif
    enddo
    
    contains
    
      real(r_kind) function intfunc(tcoord,i) result(res)
        real(r_kind)   , intent(in ) :: tcoord ! target coordinate
        integer, intent(in ) :: i
        
        res = - M(i) * ( x(i+1) - tcoord )**4 / 24. / h(i) + M(i+1) * ( tcoord - x(i) )**4 / 24. / h(i) &
              - ( y(i  ) - M(i  ) * h(i)**2 / 6. ) / h(i) * 0.5 * ( x(i+1) - tcoord )**2                &
              + ( y(i+1) - M(i+1) * h(i)**2 / 6. ) / h(i) * 0.5 * ( tcoord - x(i  ) )**2
      
      end function intfunc
  
  end subroutine spline_integration
    
  subroutine spline1(n,x,y,da,db,nt,t,ty,d1)
  !---------------------------------subroutine  comment
  !
  !  Purpose   :   三次样条之固支条件
  !    
  !-----------------------------------------------------
  !  Input  parameters  :
  !       1.   n-----插值节点个数减1，如有九个节点 则N=8
  !       2.   x ---节点自变量  为（0：N）维向量
  !       3.   y----节点因变量  （0：N）维向量
  !       4.   nt 要计算向量的维数
  !       5.   t 要计算的向量  (1:nt) 维向量
  !       6.   da  ------起点处导数 f'(x(0))
  !       7.   db -------终点处导数 f'(x(N))
  !  Output parameters  :
  !       1.   ty ---要计算的向量，（1：nt）维
  !       2.   d1 ---一阶导数，（1：nt）维
  !
  !  Common parameters  :
  !
  !----------------------------------------------------
  !  Post Script :
  !            
  !            对于要插值的向量分量，可以不必按大小排列
  !----------------------------------------------------
  
  implicit real(r_kind)(a-z)
  
  integer::n, nt
  !n为插值节点数减1，即如果有9个节点，则n=8
  !nt 要插值点的个数，及向量t,ty的维数
  
  integer::i,j,k
  
  real(r_kind)::x(0:n),y(0:n)
  
  real(r_kind)::t(nt),ty(nt)
  
  real(r_kind),optional :: d1(nt)
  
  real(r_kind)::h(0:n-1)
  
  real(r_kind)::f1(0:n-1),f2(1:n-1)
  
  real(r_kind)::u(1:n-1),lambda(1:n-1),d(0:n)
  
  real(r_kind)::M(0:n)
  
  real(r_kind)::A(0:n,0:n)
  
  do i=0,n-1
    h (i) = x(i+1) - x(i)
    f1(i) = ( y(i+1) - y(i) ) / h(i)
  end do
  
  !对固支边界条件而言   设置  d(0) 与 d(n)
  d(0) = 6./h(0  ) * (f1(0) - da     )
  d(n) = 6./h(n-1) * (db    - f1(n-1))
  
  !求得 u, lambda, d
  do i=1,n-1
    u     (i) = h(i-1) / ( h(i-1) + h(i) )
    lambda(i) = 1 - u(i)
    
    f2(i) = ( f1(i-1) - f1(i) ) / ( x(i-1) - x(i+1) )
    
    d(i) = 6. * f2(i) 
  end do
  
  !设置A矩阵值
  A = 0
  do i = 1, n-1
    a(i,i) = 2.
  end do
  
  do i = 2, n-1
    a(i,i-1) = u(i)
  end do
  
  do i = 1, n-2
    a(i,i+1) = lambda(i)
  end do
  
  !-------相比于自然条件，这里需要设置A矩阵首末行元素，其他相同
  !设置 A矩阵的首行元素
  a(0,0) = 2.
  a(0,1) = 1.
  
  !设置A矩阵元素的末行元素
  a(n,n-1) = 1.
  a(n,n  ) = 2.
  
  ! 设置右向量值
  d(1  ) = d(1  ) - u     (1  ) * M(0)
  d(n-1) = d(n-1) - lambda(n-1) * M(n)
  
  
  !call gauss(a,d,M,N+1)
  call chase(a,d,M,N+1)
      
  
  !--------以上以及求得系数
  !已经完成插值多项式的建立
  
  !------------------------------------------------
  ! 以下开始计算具体值
  do k=1,nt
    !------------
    !  对要插值向量每个分量而言，先找到其在数据中的位置
    do i=0,n-1
     if (t(k)<x(i+1)) exit 
    end do
    
    if(i==n)then
      ! 外插
      i = n - 1
    endif
    
    ty(k) = M(i) * ( x(i+1) - t(k) )**3 / 6. / h(i) +  M(i+1) * ( t(k) - x(i) )**3 / 6. / h(i)  &
          + ( y(i  ) - M(i  ) * h(i)**2 / 6. ) * ( x(i+1) - t(k) ) / h(i)                       &
          + ( y(i+1) - M(i+1) * h(i)**2 / 6. ) * ( t(k  ) - x(i) ) / h(i)
    
    if(present(d1)) d1(k) = - M(i) * (x(i+1) - t(k))**2 / (2. * h(i)) + M(i+1) * (t(k) - x(i))**2 / (2. * h(i)) &
                            + (y(i+1) - y(i)) / h(i) - (M(i+1) - M(i)) * h(i) / 6.
    
  end do

  end subroutine spline1
    
  subroutine spline3(n,x,y,nt,t,ty,d1)
    !---------------------------------subroutine  comment
    !
    !  Purpose   :   三次样条之第四类边界条件
    !    
    !-----------------------------------------------------
    !  Input  parameters  :
    !       1.   n-----插值节点个数减1，如有九个节点 则N=8
    !       2.   x ---节点自变量  为（0：N）维向量
    !       3.   y----节点因变量  （0：N）维向量
    !       4.   nt 要计算向量的维数
    !       5.   t 要计算的向量  (1:nt) 维向量
    !  Output parameters  :
    !       1.   ty ---要计算的向量，（1：nt）维
    !       2.   d1 ---一阶导数，（1：nt）维
    !
    !  Common parameters  :
    !
    !----------------------------------------------------
    !  Post Script :
    !          
    !       1.   非扭结边界条件(not-a-knot)
    !       2.   x必须从小到大排列；对于要插值的向量分量，可以不必按大小排列
    !----------------------------------------------------
    
    implicit real(r_kind)(a-z)
    
    integer::n, nt
    !n为插值节点数减1，即如果有9个节点，则n=8
    !nt 要插值点的个数，及向量t,ty的维数
    
    integer::i,j,k
    
    real(r_kind)::x(0:n),y(0:n)
    
    real(r_kind)::t(nt),ty(nt)
    
    real(r_kind),optional :: d1(nt)
    
    real(r_kind)::h(0:n-1)
    
    real(r_kind)::f1(0:n-1),f2(1:n-1)
    
    real(r_kind)::u(1:n-1),lambda(1:n-1),d(1:n-1)
    
    real(r_kind)::M(0:n),v(1:n-1)
    
    real(r_kind)::A(1:n-1,1:n-1)
    
    do i=0,n-1
      h (i) = x(i+1) - x(i)
      f1(i) = (y(i+1) - y(i)) / h(i)
    end do
    
    
    !求得 u, lambda, d
    do i=1,n-1
      u     (i) = h(i-1) / ( h(i-1) + h(i) )
      lambda(i) = 1 - u(i)
      
      f2(i) = ( f1(i-1) - f1(i) ) / ( x(i-1) - x(i+1) )
      
      d(i) = 6. * f2(i) 
    end do
    
    !设置A矩阵值
    A = 0
    do i = 2, n-2
      a(i,i) = 2.
    end do
    
    do i = 2,n-2
      a(i,i-1) = u(i)
    end do
    
    do i = 2, n-2
      a(i,i+1) = lambda(i)
    end do
    
    a(1  ,1  ) = 2. / u(1) + (h(0) + h(1)) / h(1)
    a(1  ,2  ) = lambda(1) / u(1) - h(0) / h(1)
    a(n-1,n-2) = u(n-1) / lambda(n-1) - h(n-1) / h(n-2)
    a(n-1,n-1) = 2. / lambda(n-1) + (h(n-2) + h(n-1)) / h(n-2)
    
    ! 设置右向量值
    d(1  ) = d(1  ) / u     (1  )
    d(n-1) = d(n-1) / lambda(n-1)
    
    call chase(a,d,v,N-1)
    
    do i=1,n-1
      M(i) = v(i)
    end do
    
    M(0) = (h(0  ) + h(1  )) / h(1  ) * M(1  ) - h(0  ) / h(1  ) * M(2  )
    M(n) = (h(n-2) + h(n-1)) / h(n-2) * M(n-1) - h(n-1) / h(n-2) * M(n-2)
    
    !--------以上以及求得系数
    !已经完成插值多项式的建立
    
    !------------------------------------------------
    ! 以下开始计算具体值
    do k=1,nt
      !------------
      !  对要插值向量每个分量而言，先找到其在数据中的位置
      do i=0,n-1
        if (t(k)<x(i+1)) exit
      end do
      
      if(i==n)then
        ! 外插
        i = n - 1
      endif
      
      ty(k) = M(i) * ( x(i+1) - t(k) )**3 / 6. / h(i) + M(i+1) * ( t(k) - x(i) )**3 / 6. / h(i) &
            + ( y(i  ) - M(i  ) * h(i)**2 / 6. ) * ( x(i+1) - t(k) ) / h(i)                     &
            + ( y(i+1) - M(i+1) * h(i)**2 / 6. ) * ( t(k  ) - x(i) ) / h(i)
      
      if(present(d1)) d1(k) = - M(i) * (x(i+1) - t(k))**2 / (2. * h(i)) + M(i+1) * (t(k) - x(i))**2 / (2. * h(i)) &
                              + (y(i+1) - y(i)) / h(i) - (M(i+1) - M(i)) * h(i) / 6.
    end do
  
  end subroutine spline3
  
  subroutine spline4(n,x,y,nt,t,ty,d1)
    !---------------------------------subroutine  comment
    !
    !  Purpose   :   三次样条之第二类边界条件，端点二阶导数不变
    !  END CONDITION FOR HEIGHT AND TEMPERATURE (LAPSE RATE IS CONSTANT)
    !  SM(1)=SM(2) ; SM(IMAX-1)=SM(IMAX)
    !-----------------------------------------------------
    !  Input  parameters  :
    !       1.   n-----插值节点个数减1，如有九个节点 则N=8
    !       2.   x ---节点自变量  为（0：N）维向量
    !       3.   y----节点因变量  （0：N）维向量
    !       4.   nt 要计算向量的维数
    !       5.   t 要计算的向量  (1:nt) 维向量
    !  Output parameters  :
    !       1.   ty ---要计算的向量，（1：nt）维
    !       2.   d1 ---一阶导数
    !
    !  Common parameters  :
    !
    !----------------------------------------------------
    !  Post Script :
    !          
    !       1.   自然边界条件是第一类边界条件的特殊情况
    !            本程序可以直接处理第一类边界条件
    !       2.     
    !            对于要插值的向量分量，可以不必按大小排列
    !----------------------------------------------------
    
    implicit real(r_kind)(a-z)
    
    integer::n, nt
    !n为插值节点数减1，即如果有9个节点，则n=8
    !nt 要插值点的个数，及向量t,ty的维数
    
    integer::i,j,k
    
    real(r_kind)::x(0:n),y(0:n)
    
    real(r_kind)::t(nt),ty(nt)
    
    real(r_kind),optional :: d1(nt)
    
    real(r_kind)::h(0:n-1)
    
    real(r_kind)::f1(0:n-1),f2(1:n-1)
    
    real(r_kind)::u(1:n-1),lambda(1:n-1),d(1:n-1)
    
    real(r_kind)::M(0:n),v(1:n-1)
    
    real(r_kind)::A(1:n-1,1:n-1)
    
    do i=0,n-1
      h (i) = x(i+1) - x(i)
      f1(i) = (y(i+1) - y(i)) / h(i)
    end do
    
    
    !求得 u, lambda, d
    do i=1,n-1
      u     (i) = h(i-1) / ( h(i-1) + h(i) )
      lambda(i) = 1 - u(i)
      
      f2(i) = ( f1(i-1) - f1(i) ) / ( x(i-1) - x(i+1) )
      
      d(i) = 6. * f2(i) 
    end do
    
    !设置A矩阵值
    A = 0.
    
    A(1  ,1  ) = u(1) + 2.
    A(n-1,n-1) = 2. + lambda(n-1)
    do i = 2, n-2
      a(i,i) = 2.
    end do
    
    do i = 2,n-1
      a(i,i-1) = u(i)
    end do
    
    do i = 1, n-2
      a(i,i+1) = lambda(i)
    end do
    
    call chase(a,d,v,N-1)
    
    do i=1,n-1
      M(i) = v(i)
    end do
    M(0) = M(1)
    M(n) = M(n-1)
    
    !--------以上以及求得系数
    !已经完成插值多项式的建立
    
    !------------------------------------------------
    ! 以下开始计算具体值
    do k=1,nt
      !------------
      !  对要插值向量每个分量而言，先找到其在数据中的位置
      do i=0,n-1
        if (t(k)<x(i+1)) exit
      end do
      
      if(i==n)then
        ! 外插
        i = n - 1
      endif
      
      ty(k) = M(i) * ( x(i+1) - t(k) )**3 / 6. / h(i) + M(i+1) * ( t(k) - x(i) )**3 / 6. / h(i) &
            + ( y(i  ) - M(i  ) * h(i)**2 / 6. ) * ( x(i+1) - t(k) ) / h(i)                     &
            + ( y(i+1) - M(i+1) * h(i)**2 / 6. ) * ( t(k  ) - x(i) ) / h(i)
      
      if(present(d1)) d1(k) = - M(i) * (x(i+1) - t(k))**2 / (2. * h(i)) + M(i+1) * (t(k) - x(i))**2 / (2. * h(i)) &
                              + (y(i+1) - y(i)) / h(i) - (M(i+1) - M(i)) * h(i) / 6.
    end do
    
  end subroutine spline4
  
  subroutine chase(A,f,x,N)
  !---------------------------------subroutine  comment
  !  Version   :  V1.0    
  !  Coded by  :  syz 
  !  Date      :  2010-4-9
  !-----------------------------------------------------
  !  Purpose   :  追赶法计算三对角方程组
  !              Ax=f
  !-----------------------------------------------------
  !  Input  parameters  :
  !       1.  A系数矩阵
  !       2. f 右向量
  !  Output parameters  :
  !       1.  x方程的解
  !       2.  N维数
  !  Common parameters  :
  !
  !----------------------------------------------------
  !  Post Script :
  !       1.   注意：该方法仅适用于三对角方程组
  !       2.
  !---------------------------------------------------
   
   implicit real(r_kind)(a-z)
   integer::N
   real(r_kind)::A(N,N),f(N),x(N)
   real(r_kind)::L(2:N),u(N),d(1:N-1)
   
   real(r_kind)::c(1:N-1),b(N),e(2:N)
  
   integer::i
   
   real(r_kind)::y(N)
  
  !---------------把A矩阵复制给向量e,b,c 
  do i=1,N
    b(i)=a(i,i)
  end do 
  
  do i=1,N-1
    c(i)=a(i,i+1)
  end do
  
  do i=2,N
    e(i)=a(i,i-1)
  end do
  !------------------------
  
  do i=1,N-1
    d(i)=c(i)
  end do 
  
  u(1)=b(1)
  
  do i=2,N
    L(i)=e(i)/u(i-1)
    u(i)=b(i)-L(i)*c(i-1)
  end do
  
  !------开始回带,求得y
  
  y(1)=f(1)
  do i=2,N
    y(i)=f(i)-L(i)*y(i-1)
  end do
  
  !-----开始回带，求得x
  
  x(n)=y(n)/u(n)
  
  do i=n-1,1,-1
    x(i)=(y(i)-c(i)*x(i+1))/u(i)
  end do
  
  end subroutine chase
    
  function pow(x,y)
    real(r_kind) :: pow, x, y
    
    pow = x**y
  end function pow
  
  subroutine lapack_dgbsv(A,b,x,N,kl,ku,idx,jdx)
    integer(i_kind), intent(in ) :: N  !The number of linear equations, i.e., the order of the matrix A.  N >= 0.
    integer(i_kind), intent(in ) :: kl ! The number of subdiagonals within the band of A. KL >= 0.
    integer(i_kind), intent(in ) :: ku ! The number of superdiagonals within the band of A. KU >= 0.
    real   (r_kind), intent(in ) :: A(N,N)
    real   (r_kind), intent(in ) :: b(N)
    real   (r_kind), intent(out) :: x(N)
    integer(i_kind), intent(in ), optional :: idx,jdx
    
    real   (r8) :: x_r8(N)
    real   (r8) :: AB(2 * kl + ku + 1,N) ! AB(LDA,N) ! Matrix A in band storage
    integer(i_kind) :: ipivot(N) !The pivot indices that define the permutation matrix P; row i of the matrix was interchanged with row IPIV(i).
    
    integer(i_kind) :: nRHS = 1
    integer(i_kind) :: LDA ! The leading dimension of the array AB.  LDAB >= 2*KL+KU+1.
    integer(i_kind) :: LDB ! The leading dimension of the array B.  LDB >= max(1,N).
    integer(i_kind) :: info
    
    integer(i_kind) :: i,j
    
    LDA = 2 * kl + ku + 1
    LDB = N
    
    x_r8 = b
    
    AB = 0
    do j = 1 , n
      do i = max( 1 , j - ku ) , min( n , j + kl )
        AB( kl + ku + 1 + i - j , j ) = A( i , j )
      end do
    end do
    
    call dgbsv( N, kl, ku ,NRHS, AB, LDA, ipivot, x_r8, LDB, info )
    if(info/=0)then
      print*,'Error in lapack_dgbsv with info=',info
      if(present(idx).and.present(jdx))print*,' at idx,jdx= ',idx,jdx
      ! = 0:  successful exit
      ! < 0:  if INFO  = -i, the i-th argument had an illegal value
      ! > 0:  if INFO  = i, U(i,i) is exactly zero.  The factorization
      !has been completed, but the factor U is exactly
      !singular, and the solution has not been computed.
      stop
    endif
    
    x = x_r8
    
  end subroutine lapack_dgbsv
    
  ! flip array up and down
  subroutine flipud(a)
    real(r_kind),intent(inout) :: a(:)
    real(r_kind)               :: b(lbound(a,1):ubound(a,1))
    integer                    :: i
    
    b = a
    
    do i = lbound(a,1), ubound(a,1)
      a(i) = b(ubound(a,1) - i + 1)
    enddo
  
  end subroutine flipud
    
end module math_mod
    
